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Gradients  of  the  total  energy  with  respect  to  nuclear  coordinates  within  the 

linear  combination  of  Gaussian-type  orbitals  (LCGTO)  approach  to  local-density- 

functional  theory  are  discussed.  We  explicitly  include  the  effects  of  the  fitting 

procedures  for  both  the  direct- Coulomb  and  the  exchange-correlation  energies  in 

the  evaluation  of  the  energy  gradient  expression. 
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I.  INTRODUCTION 

Accurate  local  density  functional  (LDF)  gradients  require  corrections  for  the  fact 
that  the  orbital  basis  set  is  finite1-8  and  that  the  local  one-electron  potential  is 
treated  using  a  finite  basis  set5-8.  The  linear  combination  of  Gaussian-type  orbitals 
(LCGTO)  local  density  functional  (LDF)  approach  to  molecules9,  polymers10,  and 
slabs11  is  well  established.  For  these  lower-dimensional  systems  that  are  not  periodic 
in  three  dimensions,  Fourier-transform  methods  cannot  be  used  to  readily  compute 
gradients  of  the  total  energy  (for  a  fixed  number  of  plane  waves),  unless  one  uses 
supercell  calculations  (which  are  singular  in  the  separate  fragment  limit  for  a  fixed 
number  of  plane  waves).  If  an  LDF  such  as  Xa  lends  itself  to  an  analytical  LCGTO 
treatment — i.e.  three-dimensional  numerical  integrations  or  fits  can  be  avoided— 
then  accurate  gradients  could  be  computed  given  a  converged  self-consistent-field 
(SCF)  LDF  LCGTO  total  energy5.  In  order  to  be  able  to  use  any  LDF,  a  number 
of  workers  have  recently  implemented  methods  for  computing  the  gradients  of  the 
molecular  total  energy  that  take  into  account  the  fact  that  the  orbital  basis  set  is 
finite  and  that  the  electronic  Coulomb  potential,  the  largest  part  of  the  one-electron 
potential,  is  treated  using  a  finite  expansion6-7.  For  a  different  purpose  Averill  and 
Painter12  have  considered  fitting  the  density  using  the  LCGTO  LDF  total  energy 
and  computed  exact  gradients  including  the  effects  of  the  density  fit,  assuming 
an  exact  three-dimensional  numerical  integration  of  the  exchange-correlation  (xc) 
forces.  We  consider  herein  the  effects  of  approximating  the  xc  energy  density  using  a 
necessarily  inexact  fit  to  an  LCGTO  sum  and  the  use  of  three-dimensional  numeric?-! 
integration  techniques  in  the  fitting  scheme  on  the  evaluation  of  a  gradient  with 
respect  to  nuclear  coordinates  of  the  LCGTO-LDF  total  energy  expression. 

II.  PRELIMINARIES 
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The  nonrelativistic  total  energy  and  one-electron  molecular  orbitals,  in  the  electric 
field  caused  by  nuclei  of  charges  Zk  at  positions  R,,  are  determined  by  variational 
minimization  of  the  LDF  total-energy  expression, 


E=  £ 


zkzk 


,  IR* 


-  ^  £  Jd3r  m  <f>i(r)V2<f>i(r) 


k>k<  _  1Vfc' 

P(r )  \  ,  [p\p\ 


~Y,Zk 

k 


+ 


+  (p(r)Uxe(p(r))) , 


(1) 


where  both  (/( i))  and  (/)  are  compact  notation  for  f  d3rf( r),  and  where  [/1I/2]  is 
compact  notation  for  /  d3r  1  /  d3r2/i(ri)/2(r2)/ri2.  In  this  expression  the  4 >«  are  the 
one-electron  molecular  orbitals  and  p( r), 

P(r)  =  £n,<^(r)&(r),  (2) 

t 

is  the  charge  density  when  the  zth  molecular  orbital  is  occupied  by  n,  electrons. 
(This  work  will  not  consider  the  case  where  there  is  net  spin  polarization  of  the 
electrons,  the  generalization  to  such  a  case  is  straightforward.)  The  one-electron 
molecular  orbitals  are  expressed  in  LCGTO  fashion, 


<Mr)  =  £  c#Xj(r), 


(3) 


and  satisfy  the  LDF  one-electron  secular  equation, 


j  ' 


iV72 
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v2-£ 
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+  Ki(r)  +  Vxc(p(  r)) 


|r  —  Rfc| 

Relating  the  total  energy  and  the  one-electron  equations  are  the  electrostatic  po¬ 
tential  due  to  all  electrons, 

./  P(T’ ) 


V, 


(5) 
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and  the  local  one-electron  xc  potential, 


Vxc(r) 


dp{  r) 


Uxc{p{r))  +  />(r) 


dUJA  r)) 

dp(r) 


(6) 


This  expression  can  be  inverted  to  give  the  density  derivative  of  the  xc  energy 


operator, 


dUxc  Vxc  -  Uxe 

=  (7) 

Initially,  we  follow  the  suggestion  of  Sambe  and  Felton13  and  approximate  the 
charge  density  in  Fej,  the  xc  energy  operator,  Uxc,  and  the  xc  potential,  VTC ,  using 
LCGTO  expansions, 


P( r)  ~  P(r)  =  YlfmFm(r) 

m 

Uxc(r)  «  Uxc( r)  =  UmGm(r) 

m 

Vxc(r)^Vxc(r)  =  ^mGm(r),  (8) 


centered  primarily  on  the  nuclei,  but  in  which  bond-centered  s-type  functions  can  be 
used  to  eliminate  some  high-angular-momentum  nuclear-centered  functions.  Using 
these  fits,  the  total  energy  expression9, 


E  _  y-  ZkZk' 

k>k>  1^*  — 

p(r) 


-Zzk 


lr-Rk 


k  \  i* 
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(9) 


is  stationary  with  respect  to  variations  of  the  molecular  orbital  coefficients  Cij  and 
the  charge  density  fitting  coefficients,  /m,  if  and  only  if  the  xc  fit  is  exact.  This 
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requires  the  fitting  basis  set,  Gm( r),  be  complete.  In  that  case  the  gradient  of  the 
total  energy  with  respect  to  the  Cartesian  nuclear  coordinates  is  relatively  easy  to 
evaluate6-8, 


ZkZk,(Xk  -  xk.) 


^k^k' \Xk  —  Ak>)  _  7  f  ,3  PV1  A •*-  ~  ^k) 

fek  |R.fc  -  Ri'l3  Y  kJ  |Rfc-r|3 

+  ?/m|>r‘alt|Fm|  +  ?v”(xrM*Gm))  ~a‘’{x'dx;)+c'c' 

+  Y.iM^r\-Y.UUFn\d^rl 

m  OXk  mn  dXk 


,3  p{r)(x  -  Xk) 

7  r - 

\Rk  -  r|3 


where  c.c.  means  take  the  complex  conjugate  of  the  previous  expression  within  the 
corresponding  set  of  parentheses.  In  this  expression  Pij,  is  the  matrix  element  for 
orbitals  i  and  j  of  the  one-electron  density  matrix, 

Pij  —  52  np  52  cipcjpi  (if) 

p  0 

and  Q{j  is  the  one-electron  eigenvalue  weighted  density  matrix  element 

^>j  =  52  £pnp  52  ctpc:p-  (12) 

P  ij 

The  first  two  terms  of  Eq.  10  represent  the  standard  terms  for  the  Hellmann- 
Feynman  force,  the  next  two  terms  represent  corrections  to  Hellmann- Feynman 
arising  from  the  dependence  of  the  orbital  basis  on  nuclear  coordinates1,  and  the 
last  two  terms  represent  corrections  arising  from  fitting  the  charge  density  (Eq.  8). 
Expressions  for  third  derivatives  of  the  total  energy  have  been  obtained  assuming 
that  the  xc  fitting  basis  is  complete  and  that  the  fits  are  exact8. 

The  subject  of  this  work  is  the  correction  term  that  arises  from  approximating  Uxc 
and  Vxc.  This  term  has  been  evaluated5  for  the  special  case  of  Xa,  where  both  Uxc 
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and  Vxc  are  porportional  to  p1/3,  which  can  be  treated  using  analytical  integration. 
For  all  LDF’s  this  correction  term  diminishes  with  increasing  number  of  points  and 
number  of  fitting  functions  used  to  treat  Uxc  and  I4C.6,14 


IIT  INCOMPLETE  EXCHANGE-CORRELATION  FITTING  BASES 

In  the  case  of  an  incomplete  xc  fitting  basis  set,  Eq.  10  needs  to  be  corrected  by 
adding  a  term  to  the  right-hand  side  of  the  equation, 


3EL 

dX, 


d[um  (pGm)]  dPjj 

2^  QY  \XiXj(*m) 


ijm 


dX, 


tjm 


axk 


(13) 


The  first  term  is  the  exact  differentiation  of  the  xc  energy  of  Eq.  10,  the  second  term 
comes  from  using  the  one-electron  equations,  Eq.  4,  to  eliminate  all  terms  involving 
derivatives  of  the  density  matrix  elements  from  Eq.  10,  and  the  final  term  is  the 
only  apparent  xc  term  present  in  Eq.  10. 

Because  Eq.  13  is  a  density  functional  expression,  the  terms  involving  vm  can  be 
more  compactly  expressed, 


d_El  ^3[u„(pG„>]  y-  /_dp_  \ 

dXt  dXt  ^  m\dX„  "/’ 

and  the  whole  expression  expanded, 


dEL 

dXk 


+ 


( dxk{Uxc  Vxc) ) 


* 


(14) 


(15) 


where  the  approximate  quantities  defined  in  Eq.  8  have  been  substituted  for  the 
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appropriate  expansion.  The  least-squares-fit  (LSF)  equations  for  the  xc  energy 
fitting  coefficients,  um,  must  be  differentiated  to  evaluate  this  expression. 


IV.  NUMERICAL  INTEGRATION 

In  what  follows  we  will  assume  that  we  have  at  our  disposal  a  numerical  integration 
scheme  involving  a  set  of  points,  r,  and  an  associated  set  of  volume  elements  At,, 

{/(**)}=  £  Ar«  /(r.),  (16) 

i 

with  the  property  that  we  can  increase  the  number  of  points  to  make  the  numerical 
integral  {/(r)}  arbitrarily  close  in  value  to  the  analytic  integral  (/( r)).  We  further 
assume  that  numerical  integration  is  used  only  to  fit  Uxc  and  Vxc ,  via  a  weighted 
LSF  scheme, 

£vm{WvGmGn}  =  {WvVxcGn},  (17) 

m 

£  um  { WuGmGn]  =  {WuUxcGn} ,  (18) 


where  the  weight  functions  Wv  and  Wu  axe  arbitrary  positive  (or  negative)  semidef- 
inite  functions  of  position.  Eq.  18  can  be  inverted  to  give  the  fitting  coefficients, 

um  =  £  {WuGG}^  {WuUxcGn}  ,  (19) 

n 

where  the  symmetric  matrix  GG  is  the  tensor  product  of  the  vector  G  with  itself. 

We  now  assume  that  we  have  some  method  for  choosing  the  points  and  volume 
elements  of  Eq.  16  so  that  we  can  systematically  improve  the  accuracy  of  this 
numerical  integration,  albeit  at  an  ever  increasing  computational  price.  Therefore, 
we  shall  ignore  how  the  weights  and  points  of  this  numerical  integration  vary  as  any 
Cartesian  nuclear  coordinate  changes. 
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Differentiating  Eq.  18  gives  the  derivatives  of  the  xc  energy  fitting  coefficients, 
+  | §£(*«  -  fUGn}  +  lwu(Uxc  -  Uxc)^X\  .  (2( 


The  derivatives  of  the  xc  energy  fitting  coefficients  in  Eq.  15  appear  multiplied  by 
the  overlap  integral  of  the  density  and  the  xc  basis  functions  suggesting  a  new  set 
of  variables, 

(™  =  E(,l/uGG  U(pG„>,  (21) 


T(r)  =  £tmG„(r). 


Using  these  variables,  the  summation, 


+  {  -  (7«)r|  +  ■£  u  |  wa(uce  - , 


can  be  used  to  eliminate  the  derivatives  of  the  xc  energy  density  fitting  coefficients 


from  Eq.  15, 


If  =  {If v~  ~  u~v)  +  {"W-  -  ^if} 

<24> 
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where  Eq.  7  has  been  used  to  evaluate  the  derivative  of  Uxc. 

Eqs.  21  and  22  can  be  rewritten, 

{WUTG}  =  (pG).  (25) 

In  other  words,  for  each  weight  function  the  corresponding  T  exactly  transfoims 
a  numerical  integral  into  an  analytic  integral  on  the  set  of  functions  exactly  fit  by 
the  basis  G.  For  other  functions,  the  transformation  is  into  the  analytic  integral  of 
the  fitted  quantity, 

{WUTF}  =  (, pF ).  (26) 

In  particular,  this  hold  for  the  xc  energy  operator, 

{WUTUXC}  =  (pUxc),  (27) 

and 

{WUTVXC}  =  (pVxc),  (28) 

provided  Wv  is  chosen  to  be  Wu.  These  relationships  are  independent  of  the  accuracy 
of  the  numerical  integration  scheme  and  requires  only  that  the  matrix  {W^GG}  be 
invertible.  If  the  basis  G  is  complete,  then  Eq.  23  vanishes  identically  independent  of 
the  quality  of  the  numerical  integration.  Therefore,  for  large  xc  basis  sets  improving 
the  numerical  integration  scheme  might  not  improve  the  accuracy  of  the  gradients 
if  they  are  approximately  computed  using  Eq.  10. 

V.  VARIATIONAL  XC  FITTING 

It  appears  impossible  to  avoid  resolving  the  one-electron  equations  to  determine 
derivatives  of  p  with  respect  to  the  nuclear  coordinates  that  are  necessary  to  eval¬ 
uate  Eq.  23  precisely  when  the  fitting  basis  set  is  incomplete.  An  additional  SCF 
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processes,  analogous  to  the  coupled-perturbed-Hartree-Fock  calculation15,  can  be 
avoided  if  the  one-electron  equations  are  true  variations  of  Eq.  9.  To  satisfy  this 
requirement,  a  variational  xc  potential  must  be  defined,  perhaps  most  directly, 
by  considering  the  derivative  with  respect  to  an  arbitrary  one-electron  occupation 
number, 


(^■) 


d(pUxc) 

dr>i 

=  (GU^)  +  E 

m  unt 


(29) 


The  remaining  derivatives  can  ’-e  obtained  by  differentiating  Eq.  18  again.  This 
equation  contains  half  as  many  terms  as  Eq.  24  because  G  is  now  assumed  constant, 


(<t>:vxr<t>i)  =  +  \  4>:nuxc  -  u 


dp 


(30) 


In  this  equation  the  weight  function  does  not  contain  a  subscript,  because  the  xc 
potential  is  not  going  to  be  fit.  Instead  the  matrix  elements  of  the  xc  potential, 

{x'XPxi)  =  (x- u,'Xj) 

+  jx'T  (([/,.  -  +  (Vxc  -  Ui.Jj)  *,}  .  (31) 

were  obtained  directly.  In  this  expression  W  is  taken  to  be  a  function  of  the  density 
alone. 

If  the  weight  function  is  chosen  to  be  unity,  then  summing  this  expression  over 
all  occupied  orbitals  and  using  Eqs.  27  and  28  gives  precisely  {pVxc),  as  expected. 
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The  weight  function16, 


w  =  JT’  (32> 

that  is  the  optimal  density  functional  expression  for  computing  the  total  energy  in 
the  nonvariational  LCGTO  method9-11  has  derivative 


dW  =  2  Uxc  -  VIC 

dp  Ul 


(33) 


The  choice  of  weighting  function,  W,  is  less  important  now  because  the  function  is 
included  in  the  SCF  process. 

Having  defined  a  variational  xc  potential,  the  derivative  of  the  fitted  xc  energy  is 
simplified, 


where  the  parentheses  and  subscript  on  the  derivative  of  the  xc  energy  coefficient 
indicates  that  p  must  be  held  constant  during  this  partial  differentiation.  The 
first  term  comes  from  Eq.  29,  and  the  remaining  two  terms  come  from  changes  in 
the  fitting  basis  set,  which  is  independent  of  the  density.  Using  the  variational  xc 
potential,  Eq.  31,  in  the  one-electron  secular  equation,  Eq.  4,  eliminates  all  terms 
containing  a  derivative  of  the  density  (and  density  functional  weight  function)  from 
the  gradient  correction  term,  Eq.  24.  Thus  the  precise  derivative  of  the  total  energy 
using  a  variational  xc  potential  is 
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where 


:3at 


r  /  _  (/w  w 

+  I x’T  f (Uxc  -  Uxc)—  +  (Vxe  -  Uxc)— 


dXL 

dXk 


(36) 


The  price  that  is  paid  for  an  exact  derivative  is  the  coupling  of  the  eigenvalue 
problem  and  numerical  fitting/integration. 

VI.  CONCLUSION 

The  error  terms  in  existing  gradient  LCGTO  LDF  programs  have  been  analyzed. 
The  derivatives  of  the  fitting  coefficients  with  respect  to  nuclear  coordinate  have 
been  eliminated  from  Eq.  23,  however  the  derivative  of  the  density  cannot  be  elim¬ 
inated  without  adding  an  additional  SCF  proceedure.  A  better  approach  is  to  use 
Vxvcar  in  the  one-electron  equations.  This  correction  replace  two  quasi-independent 
fitting  proceedures,  Eq.  17  and  Eq.  18,  with  one,  the  latter.  Thus  the  xc  fitting 
basis  set  can  be  optimized  to  fit  only  the  xc  energy  operator.  Furthermore,  this 
correction  leads  to  directly  computable  gradients,  Eq.  35.  This  derivative  is  exact 
under  the  assumption  that  the  grid  of  points  necessary  to  fit  the  xc  energy  density 
do  not  move  with  the  nuclei.  That  is,  of  course,  not  true  unless  the  grid  leads  to 
accurate  numerical  integration.  Eq.  26  strongly  suggests  that  this  remaining  prob¬ 
lem  will  be  minimized.  In  any  event,  existing  grids17  make  this  problem  negligible, 
at  least  temporarily. 
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